sampleDists <- dist(t(assay(rld)))
plot(hclust(sampleDists))annot = dplyr::select(experimental_metadata, condition)
row.names(annot) = experimental_metadata$sample_id
rld %>%
assay() %>%
cor() %>%
pheatmap(color = colorRampPalette(rev(brewer.pal(n = 7, name = "RdBu")))(100),
annotation = annot,
annotation_colors = list(condition = c(Control = "#BF3100", OKSM = "#E9B44C",
Senolytic = "#1B998B", Senolytic_OKSM = "#5D576B")),
cluster_rows = TRUE,
cluster_cols = T,
cellwidth = 13,
cellheight = 13)ntop = 500
rv <- rowVars(assay(rld))
select <- order(rv, decreasing = TRUE)[seq_len(min(ntop, length(rv)))]
pca = prcomp(t(assay(rld)[select,]))
percentVar <- pca$sdev^2/sum(pca$sdev^2)
pca_data <- plotPCA(rld, intgroup = c("condition", "replicate"), returnData=TRUE)
percentVar <- round(100 * attr(pca_data, "percentVar"), digits=2)
ggplot(pca_data, aes(PC1, PC2, color=condition, shape=replicate)) + geom_point(size=3) +
scale_x_continuous(paste0("PC1: ",percentVar[1],"% variance")) +
scale_y_continuous(paste0("PC2: ",percentVar[2],"% variance")) +
coord_fixed() + theme_classic() + geom_text(data = pca_data, aes(PC1,PC2, label = name), hjust = 1.2)From the hierarchical clustering, correlation heatmap, and PCA, there seems to be some variation that we’re not accounting for. Here, we use SVA to account for the unknown variation. Since the number of variables is unknown, we call the sva function without the n.sv argument, allowing the algorithm to estimate the number of factors.
# There are batch effects for this dataset. However, running ~ batch + condition gives us an error that says the model matrix is not full rank. So we follow this guide https://master.bioconductor.org/packages/release/workflows/vignettes/rnaseqGene/inst/doc/rnaseqGene.html#removing-hidden-batch-effects
library(sva)
dat <- counts(dds, normalized = TRUE)
idx <- rowMeans(dat) > 1
dat <- dat[idx, ]
mod <- model.matrix(~ condition, colData(rld))
mod0 <- model.matrix(~ 1, colData(rld))
svseq <- svaseq(dat, mod, mod0)## Number of significant surrogate variables is: 3
## Iteration (out of 5 ):1 2 3 4 5
According to SVA, there are 3 significant surrogate variables.
par(mfrow = c(3, 1), mar = c(3,5,3,1))
for (i in 1:3) {
stripchart(svseq$sv[, i] ~ dds$condition, vertical = TRUE, main = paste0("SV", i))
abline(h = 0)
}After accounting for the 3 significant surrogate variables, this is what the data looks like:
sampleDists <- dist(t(assay(rld)))
plot(hclust(sampleDists))rld %>%
assay() %>%
cor() %>%
pheatmap(color = colorRampPalette(rev(brewer.pal(n = 7, name = "RdBu")))(100),
annotation = annot,
annotation_colors = list(condition = c(Control = "#BF3100", OKSM = "#E9B44C",
Senolytic = "#1B998B", Senolytic_OKSM = "#5D576B")),
cluster_rows = TRUE,
cluster_cols = T,
cellwidth = 13,
cellheight = 13) The correlation heatmap shows that the sample AT_Sen_1 could be an outlier.
ntop = 500
rv <- rowVars(assay(rld))
select <- order(rv, decreasing = TRUE)[seq_len(min(ntop, length(rv)))]
pca = prcomp(t(assay(rld)[select,]))
percentVar <- pca$sdev^2/sum(pca$sdev^2)
pca_data <- plotPCA(rld, intgroup = c("condition", "replicate"), returnData=TRUE)
percentVar <- round(100 * attr(pca_data, "percentVar"), digits=2)
ggplot(pca_data, aes(PC1, PC2, color=condition, shape=replicate)) + geom_point(size=3) +
scale_x_continuous(paste0("PC1: ",percentVar[1],"% variance")) +
scale_y_continuous(paste0("PC2: ",percentVar[2],"% variance")) +
coord_fixed() + theme_classic() + geom_text(data = pca_data, aes(PC1,PC2, label = name), hjust = 1.2)Taking the three plots into account, AT_Sen_1 is removed.
We re-run sva on the new dataset not including AT_Sen_1:
# There are batch effects for this dataset. However, running ~ batch + condition gives us an error that says the model matrix is not full rank. So we follow this guide https://master.bioconductor.org/packages/release/workflows/vignettes/rnaseqGene/inst/doc/rnaseqGene.html#removing-hidden-batch-effects
library(sva)
dat <- counts(dds, normalized = TRUE)
idx <- rowMeans(dat) > 1
dat <- dat[idx, ]
mod <- model.matrix(~ condition, colData(rld))
mod0 <- model.matrix(~ 1, colData(rld))
svseq <- svaseq(dat, mod, mod0)## Number of significant surrogate variables is: 2
## Iteration (out of 5 ):1 2 3 4 5
According to SVA, there are 2 significant surrogate variables.
par(mfrow = c(2, 1), mar = c(3,5,3,1))
for (i in 1:2) {
stripchart(svseq$sv[, i] ~ dds$condition, vertical = TRUE, main = paste0("SV", i))
abline(h = 0)
}After accounting for the 2 significant surrogate variables, this is what the data looks like:
sampleDists <- dist(t(assay(rld)))
plot(hclust(sampleDists))rld %>%
assay() %>%
cor() %>%
pheatmap(color = colorRampPalette(rev(brewer.pal(n = 7, name = "RdBu")))(100),
annotation = annot,
annotation_colors = list(condition = c(Control = "#BF3100", OKSM = "#E9B44C",
Senolytic = "#1B998B", Senolytic_OKSM = "#5D576B")),
cluster_rows = TRUE,
cluster_cols = T,
cellwidth = 13,
cellheight = 13)ntop = 500
rv <- rowVars(assay(rld))
select <- order(rv, decreasing = TRUE)[seq_len(min(ntop, length(rv)))]
pca = prcomp(t(assay(rld)[select,]))
percentVar <- pca$sdev^2/sum(pca$sdev^2)
pca_data <- plotPCA(rld, intgroup = c("condition", "replicate"), returnData=TRUE)
percentVar <- round(100 * attr(pca_data, "percentVar"), digits=2)
ggplot(pca_data, aes(PC1, PC2, color=condition, shape=replicate)) + geom_point(size=3) +
scale_x_continuous(paste0("PC1: ",percentVar[1],"% variance")) +
scale_y_continuous(paste0("PC2: ",percentVar[2],"% variance")) +
coord_fixed() + theme_classic() + geom_text(data = pca_data, aes(PC1,PC2, label = name), hjust = 1.2)effective_lengths = matrix(0, ncol=length(experimental_metadata$sample_id), nrow=17714)
colnames(effective_lengths)= experimental_metadata$sample_id
for( i in experimental_metadata$sample_id){
effective_lengths[,i] = read.table(paste("data/ubiquitous/", i, ".genes.results",sep=""), sep="\t", header=TRUE)$effective_length
}
row.names(effective_lengths) = read.table(paste("data/ubiquitous/", i, ".genes.results",sep=""), sep="\t", header=TRUE)$gene_id
effective_lengths = rowMeans(effective_lengths[row.names(counts(ddssva)),])
ncrpk = counts(ddssva) / (effective_lengths / 1000)
ncrpk = apply(ncrpk, c(1,2), function(x){if(is.nan(x)){0}else{x}})
ncrpk = apply(ncrpk, c(1,2), function(x){if(is.infinite(x)){0}else{x}})
ncscalingfactor = colSums(ncrpk) / 1e6
nctpm = sweep(ncrpk, 2, ncscalingfactor, "/")
nctpm.melt = melt(nctpm)
ggplot(nctpm.melt, aes(x=Var2, y=value)) + geom_boxplot() + theme_classic() + theme(axis.text.x = element_text(angle = 90, colour="black", hjust = 1)) + scale_x_discrete("Sample") + scale_y_continuous("TPM")tpm.threshold = 10000
test.tpm = apply(nctpm, 1, function(x){ any(x> tpm.threshold) })
ensembl.genes[names(test.tpm[test.tpm])] %>%
as.data.frame() %>%
datatable(options = list(scrollX = TRUE))generate_de_section(dds_wald, "OKSM", "Control")## [1] "Number of significant genes (padj < 0.1 & log2FoldChange < log2(1.5)): 317"
generate_de_section(dds_wald, "Senolytic", "Control")## [1] "Number of significant genes (padj < 0.1 & log2FoldChange < log2(1.5)): 550"
generate_de_section(dds_wald, "Senolytic_OKSM", "Control")## [1] "Number of significant genes (padj < 0.1 & log2FoldChange < log2(1.5)): 284"
generate_de_section(dds_wald, "OKSM", "Senolytic")## [1] "Number of significant genes (padj < 0.1 & log2FoldChange < log2(1.5)): 460"
generate_de_section(dds_wald, "OKSM", "Senolytic_OKSM")## [1] "Number of significant genes (padj < 0.1 & log2FoldChange < log2(1.5)): 170"
generate_de_section(dds_wald, "Senolytic", "Senolytic_OKSM")## [1] "Number of significant genes (padj < 0.1 & log2FoldChange < log2(1.5)): 493"
dds_LRT = nbinomLRT(ddssva, reduced = ~SV1+SV2)
res = results(dds_LRT)
res$gene_biotype= ensembl.genes$gene_biotype[match(row.names(res), ensembl.genes$gene_id)]
res$external_gene_name= ensembl.genes$external_gene_name[match(row.names(res), ensembl.genes$gene_id)]
hist(res$pvalue)Number of significant genes (padj < 0.1):
sum(res$padj < 0.1, na.rm=T)## [1] 1282
# assay(x) to access the count data
assay(significant_rld) <- limma::removeBatchEffect(assay(significant_rld), covariates = cov)
sig_mat_rld = assay(significant_rld)
# The apply function swaps the rows to samples and columns to genes -- the standard is the other way around: samples in cols and genes in rows, hence the transpose function
zscores = t(apply(sig_mat_rld, 1, function(x){ (x - mean(x)) / sd(x) }))foo = as(zscores, "matrix")
bar = sapply(1:20, function(x){kmeans(foo, centers=x)$tot.withinss})
plot(bar, type="l")pam_clust <- generate_data(zscores, 7, "pam")
saveRDS(pam_clust, "output/ubiquitous/pam_clust.rds")
# pam_clust <- as.data.frame(pam_clust)
# pam_clust$Cluster <- factor(pam_clust$Cluster, levels = c(5,4,3,1,2,6))
# pam_clust <- pam_clust[order(pam_clust$Cluster),]
pheatmap(pam_clust[,1:(ncol(pam_clust)-1)],
color = colorRampPalette(rev(brewer.pal(n = 7, name = "RdBu")))(100),
fontsize_row = 5.5,
annotation_col = annotation,
annotation_colors = anno_colours,
cluster_rows = FALSE,
cluster_cols = FALSE)pam_clust_df <- pam_clust %>%
as.data.frame() %>%
mutate(gene_name = ensembl.genes[rownames(.),]$external_gene_name) %>%
dplyr::select(gene_name, Cluster) %>%
dplyr::rename("Gene Name" = gene_name)
datatable(pam_clust_df, options = list(scrollX = TRUE), class = "white-space: nowrap")c1 <- pam_clust_df[pam_clust_df$Cluster == 1, ] %>%
dplyr::select(-Cluster)
c2 <- pam_clust_df[pam_clust_df$Cluster == 2, ] %>%
dplyr::select(-Cluster)
c3 <- pam_clust_df[pam_clust_df$Cluster == 3, ] %>%
dplyr::select(-Cluster)
c4 <- pam_clust_df[pam_clust_df$Cluster == 4, ] %>%
dplyr::select(-Cluster)
c5 <- pam_clust_df[pam_clust_df$Cluster == 5, ] %>%
dplyr::select(-Cluster)
c6 <- pam_clust_df[pam_clust_df$Cluster == 6, ] %>%
dplyr::select(-Cluster)
c7 <- pam_clust_df[pam_clust_df$Cluster == 7, ] %>%
dplyr::select(-Cluster)
data.frame(Cluster = c("Cluster 1", "Cluster 2", "Cluster 3", "Cluster 4", "Cluster 5",
"Cluster 6", "Cluster 7", "Total"),
Number_of_genes = c(nrow(c1), nrow(c2), nrow(c3), nrow(c4), nrow(c5),
nrow(c6), nrow(c7),
sum(c(nrow(c1),nrow(c2),nrow(c3),nrow(c4),nrow(c5),
nrow(c6),nrow(c7)))))## Cluster Number_of_genes
## 1 Cluster 1 300
## 2 Cluster 2 201
## 3 Cluster 3 157
## 4 Cluster 4 173
## 5 Cluster 5 147
## 6 Cluster 6 204
## 7 Cluster 7 100
## 8 Total 1282
## Checking the stability of clusters
#fviz_nbclust(zscores, kmeans, method = "silhouette")
set.seed(2)
dd = as.dist((1 - cor(t(zscores)))/2)
pam = pam(dd, 7, diss = TRUE)
sil = silhouette(pam$clustering, dd)
plot(sil, border=NA, main = "Silhouette Plot for 7 Clusters")#listEnrichrSites()
setEnrichrSite("FlyEnrichr")
dbs <- listEnrichrDbs()
to_check <- c("GO_Biological_Process_2018", "KEGG_2019")
output_dir <- "output/ubiquitous/"eresList$GO_Biological_Process_2018 %>%
plot_enrichr("GO_Biological_Process_2018")datatable(eresList[[1]], options = list(scrollX = TRUE), class = "white-space: nowrap")write.csv(eresList$GO_Biological_Process_2018, paste0(output_dir,"c1_go.csv"), row.names = T, col.names = T)
get_wanted_terms(eresList$GO_Biological_Process_2018, output_dir, "ubi_c1_go", c("mitotic cytokinesis (GO:0000281)", "DNA replication (GO:0006260)"))eresList$KEGG_2019 %>%
plot_enrichr("KEGG_2019") datatable(eresList[[2]], options = list(scrollX = TRUE), class = "white-space: nowrap")write.csv(eresList$KEGG_2019, paste0(output_dir,"c1_kegg.csv"), row.names = T, col.names = T)
get_wanted_terms(eresList$KEGG_2019, output_dir, "ubi_c1_kegg", c("Base excision repair", "FoxO signaling pathway"))eresList$GO_Biological_Process_2018 %>%
plot_enrichr("GO_Biological_Process_2018")datatable(eresList[[1]], options = list(scrollX = TRUE), class = "white-space: nowrap")write.csv(eresList$GO_Biological_Process_2018, paste0(output_dir,"c2_go.csv"), row.names = T, col.names = T)
get_wanted_terms(eresList$GO_Biological_Process_2018, output_dir, "ubi_c2_go", c("glutathione metabolic process (GO:0006749)", "peptide metabolic process (GO:0006518)"))eresList$KEGG_2019 %>%
plot_enrichr("KEGG_2019") datatable(eresList[[2]], options = list(scrollX = TRUE), class = "white-space: nowrap")write.csv(eresList$KEGG_2019, paste0(output_dir,"c2_kegg.csv"), row.names = T, col.names = T)
get_wanted_terms(eresList$KEGG_2019, output_dir, "ubi_c2_kegg", c("Starch and sucrose metabolism","Lysosome"))eresList$GO_Biological_Process_2018 %>%
plot_enrichr("GO_Biological_Process_2018")datatable(eresList[[1]], options = list(scrollX = TRUE), class = "white-space: nowrap")write.csv(eresList$GO_Biological_Process_2018, paste0(output_dir,"c3_go.csv"), row.names = T, col.names = T)
get_wanted_terms(eresList$GO_Biological_Process_2018, output_dir, "ubi_c3_go", c("proteolysis (GO:0006508)"))eresList$KEGG_2019 %>%
plot_enrichr("KEGG_2019") datatable(eresList[[2]], options = list(scrollX = TRUE), class = "white-space: nowrap")write.csv(eresList$KEGG_2019, paste0(output_dir,"c3_kegg.csv"), row.names = T, col.names = T)
get_wanted_terms(eresList$KEGG_2019, output_dir, "ubi_c3_kegg", c("Lysosome", "Amino sugar and nucleotide sugar metabolism"))eresList$GO_Biological_Process_2018 %>%
plot_enrichr("GO_Biological_Process_2018")datatable(eresList[[1]], options = list(scrollX = TRUE), class = "white-space: nowrap")write.csv(eresList$GO_Biological_Process_2018, paste0(output_dir,"c4_go.csv"), row.names = T, col.names = T)
get_wanted_terms(eresList$GO_Biological_Process_2018, output_dir, "ubi_c4_go", c("phagocytosis (GO:0006909)", "secretion (GO:0046903)"))eresList$KEGG_2019 %>%
plot_enrichr("KEGG_2019") datatable(eresList[[2]], options = list(scrollX = TRUE), class = "white-space: nowrap")write.csv(eresList$KEGG_2019, paste0(output_dir,"c4_kegg.csv"), row.names = T, col.names = T)
get_wanted_terms(eresList$KEGG_2019, output_dir, "ubi_c4_kegg", c("Pantothenate and CoA biosynthesis", "Neuroactive ligand-receptor interaction"))eresList$GO_Biological_Process_2018 %>%
plot_enrichr("GO_Biological_Process_2018")datatable(eresList[[1]], options = list(scrollX = TRUE), class = "white-space: nowrap")write.csv(eresList$GO_Biological_Process_2018, paste0(output_dir,"c5_go.csv"), row.names = T, col.names = T)
get_wanted_terms(eresList$GO_Biological_Process_2018, output_dir, "ubi_c5_go", c("translation (GO:0006412)"))eresList$KEGG_2019 %>%
plot_enrichr("KEGG_2019") datatable(eresList[[2]], options = list(scrollX = TRUE), class = "white-space: nowrap")write.csv(eresList$KEGG_2019, paste0(output_dir,"c5_kegg.csv"), row.names = T, col.names = T)
get_wanted_terms(eresList$KEGG_2019, output_dir, "ubi_c5_kegg", c("Starch and sucrose metabolism", "Ribosome"))eresList$GO_Biological_Process_2018 %>%
plot_enrichr("GO_Biological_Process_2018")datatable(eresList[[1]], options = list(scrollX = TRUE), class = "white-space: nowrap")write.csv(eresList$GO_Biological_Process_2018, paste0(output_dir,"c6_go.csv"), row.names = T, col.names = T)
get_wanted_terms(eresList$GO_Biological_Process_2018, output_dir, "ubi_c6_go", c("copper ion homeostasis (GO:0055070)", "chloride transport (GO:0006821)"))eresList$KEGG_2019 %>%
plot_enrichr("KEGG_2019") datatable(eresList[[2]], options = list(scrollX = TRUE), class = "white-space: nowrap")write.csv(eresList$KEGG_2019, paste0(output_dir,"c6_kegg.csv"), row.names = T, col.names = T)
get_wanted_terms(eresList$KEGG_2019, output_dir, "ubi_c6_kegg", c("Glutathione metabolism", "Insect hormone biosynthesis"))eresList$GO_Biological_Process_2018 %>%
plot_enrichr("GO_Biological_Process_2018")datatable(eresList[[1]], options = list(scrollX = TRUE), class = "white-space: nowrap")write.csv(eresList$GO_Biological_Process_2018, paste0(output_dir,"c7_go.csv"), row.names = T, col.names = T)
get_wanted_terms(eresList$GO_Biological_Process_2018, output_dir, "ubi_c7_go", c("regulation of cell proliferation (GO:0042127)", "Notch signaling pathway (GO:0007219)"))eresList$KEGG_2019 %>%
plot_enrichr("KEGG_2019") datatable(eresList[[2]], options = list(scrollX = TRUE), class = "white-space: nowrap")write.csv(eresList$KEGG_2019, paste0(output_dir,"c7_kegg.csv"), row.names = T, col.names = T)
get_wanted_terms(eresList$KEGG_2019, output_dir, "ubi_c7_kegg", c("Arginine and proline metabolism"))ame_res[[1]] %>%
dplyr::select(-motif_DB, -motif_ID) %>%
datatable(options = list(scrollX = TRUE))ame_res[[2]] %>%
dplyr::select(-motif_DB, -motif_ID) %>%
datatable(options = list(scrollX = TRUE))ame_res[[3]] %>%
dplyr::select(-motif_DB, -motif_ID) %>%
datatable(options = list(scrollX = TRUE))ame_res[[4]] %>%
dplyr::select(-motif_DB, -motif_ID) %>%
datatable(options = list(scrollX = TRUE))ame_res[[5]] %>%
dplyr::select(-motif_DB, -motif_ID) %>%
datatable(options = list(scrollX = TRUE))ame_res[[6]] %>%
dplyr::select(-motif_DB, -motif_ID) %>%
datatable(options = list(scrollX = TRUE))ame_res[[7]] %>%
dplyr::select(-motif_DB, -motif_ID) %>%
datatable(options = list(scrollX = TRUE))sessionInfo()## R version 4.0.5 (2021-03-31)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Big Sur 10.16
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] stats4 parallel stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] sva_3.38.0 BiocParallel_1.24.1
## [3] mgcv_1.8-34 nlme_3.1-152
## [5] factoextra_1.0.7 enrichR_3.0
## [7] DT_0.17 cowplot_1.1.1
## [9] RColorBrewer_1.1-2 scales_1.1.1
## [11] reshape2_1.4.4 knitr_1.33
## [13] biomaRt_2.46.3 GenomicFeatures_1.42.3
## [15] AnnotationDbi_1.52.0 genefilter_1.72.1
## [17] DESeq2_1.30.1 SummarizedExperiment_1.20.0
## [19] Biobase_2.50.0 MatrixGenerics_1.2.1
## [21] matrixStats_0.58.0 BSgenome.Dmelanogaster.UCSC.dm6_1.4.1
## [23] BSgenome_1.58.0 rtracklayer_1.50.0
## [25] GenomicRanges_1.42.0 GenomeInfoDb_1.26.7
## [27] Biostrings_2.58.0 XVector_0.30.0
## [29] IRanges_2.24.1 S4Vectors_0.28.1
## [31] BiocGenerics_0.36.0 forcats_0.5.1
## [33] stringr_1.4.0 dplyr_1.0.5
## [35] purrr_0.3.4 readr_1.4.0
## [37] tidyr_1.1.3 tibble_3.1.0
## [39] tidyverse_1.3.0 EnhancedVolcano_1.8.0
## [41] ggrepel_0.9.1 ggplot2_3.3.3
## [43] pheatmap_1.0.12 cluster_2.1.1
##
## loaded via a namespace (and not attached):
## [1] readxl_1.3.1 backports_1.2.1 BiocFileCache_1.14.0
## [4] plyr_1.8.6 splines_4.0.5 crosstalk_1.1.1
## [7] digest_0.6.27 htmltools_0.5.1.1 fansi_0.4.2
## [10] magrittr_2.0.1 memoise_2.0.0 openxlsx_4.2.3
## [13] limma_3.46.0 annotate_1.68.0 modelr_0.1.8
## [16] extrafont_0.17 extrafontdb_1.0 askpass_1.1
## [19] prettyunits_1.1.1 colorspace_2.0-0 blob_1.2.1
## [22] rvest_1.0.0 rappdirs_0.3.3 haven_2.3.1
## [25] xfun_0.22 crayon_1.4.1 RCurl_1.98-1.3
## [28] jsonlite_1.7.2 survival_3.2-10 glue_1.4.2
## [31] gtable_0.3.0 zlibbioc_1.36.0 DelayedArray_0.16.3
## [34] proj4_1.0-10.1 car_3.0-10 Rttf2pt1_1.3.8
## [37] maps_3.3.0 abind_1.4-5 edgeR_3.32.1
## [40] DBI_1.1.1 rstatix_0.7.0 Rcpp_1.0.6
## [43] xtable_1.8-4 progress_1.2.2 foreign_0.8-81
## [46] bit_4.0.4 htmlwidgets_1.5.3 httr_1.4.2
## [49] ellipsis_0.3.1 farver_2.1.0 pkgconfig_2.0.3
## [52] XML_3.99-0.6 sass_0.3.1 dbplyr_2.1.1
## [55] locfit_1.5-9.4 utf8_1.2.1 labeling_0.4.2
## [58] tidyselect_1.1.0 rlang_0.4.10 munsell_0.5.0
## [61] cellranger_1.1.0 tools_4.0.5 cachem_1.0.4
## [64] cli_2.4.0 generics_0.1.0 RSQLite_2.2.6
## [67] broom_0.7.6 evaluate_0.14 fastmap_1.1.0
## [70] yaml_2.2.1 bit64_4.0.5 fs_1.5.0
## [73] zip_2.1.1 ash_1.0-15 ggrastr_0.2.3
## [76] xml2_1.3.2 compiler_4.0.5 rstudioapi_0.13
## [79] beeswarm_0.3.1 curl_4.3 ggsignif_0.6.1
## [82] reprex_2.0.0 geneplotter_1.68.0 bslib_0.2.4
## [85] stringi_1.5.3 highr_0.8 ggalt_0.4.0
## [88] lattice_0.20-41 Matrix_1.3-2 vctrs_0.3.7
## [91] pillar_1.6.0 lifecycle_1.0.0 jquerylib_0.1.3
## [94] data.table_1.14.0 bitops_1.0-6 R6_2.5.0
## [97] rio_0.5.26 KernSmooth_2.23-18 vipor_0.4.5
## [100] MASS_7.3-53.1 assertthat_0.2.1 rjson_0.2.20
## [103] openssl_1.4.3 rprojroot_2.0.2 withr_2.4.1
## [106] GenomicAlignments_1.26.0 Rsamtools_2.6.0 GenomeInfoDbData_1.2.4
## [109] hms_1.0.0 grid_4.0.5 rmarkdown_2.7
## [112] carData_3.0-4 ggpubr_0.4.0 lubridate_1.7.10
## [115] ggbeeswarm_0.6.0